Obtención de los registros meteorológicos del SENAMHI

La captura de datos fueron obtenidos de la pagina correspondiente a la descarga del SENAMHI.

Información registrada

“year”, “mes”, “dia”, “prec”

1. Integración de datos

1.1 Lectura de las estaciones

# Enumerando las estaciones
L01<-read.table("Ayabaca.txt"); L01$estacion=1
L02<-read.table("Bernal.txt"); L02$estacion=2
L03<-read.table("Chulucanas.txt"); L03$estacion=3
L04<-read.table("Chusis.txt"); L04$estacion=4
L05<-read.table("Hacienda Bigote.txt"); L05$estacion=5
L06<-read.table("Hacienda Sumaya.txt"); L06$estacion=6
L07<-read.table("Huarmaca.txt"); L07$estacion=7
L08<-read.table("La Esperanza.txt"); L08$estacion=8
L09<-read.table("Lancones.txt"); L09$estacion=9
L10<-read.table("Mallares.txt"); L10$estacion=10
L11<-read.table("Miraflores.txt"); L11$estacion=11
L12<-read.table("Morropon.txt"); L12$estacion=12
L13<-read.table("Pananga.txt"); L13$estacion=13
L14<-read.table("Porculla.txt"); L14$estacion=14
L15<-read.table("San Miguel.txt"); L15$estacion=15
L16<-read.table("San Pedro.txt"); L16$estacion=16
L17<-read.table("Sapillica.txt"); L17$estacion=17
L18<-read.table("Sausal de Culucan.txt"); L18$estacion=18
L19<-read.table("Sondorillo.txt"); L19$estacion=19
L20<-read.table("Tuluce.txt"); L20$estacion=20

1.2 Integración y Limpieza de datos

data_anio=function(anio){
  R01<-L01[(L01[,1]==anio),]
  R02<-L02[(L02[,1]==anio),]
  R03<-L03[(L03[,1]==anio),]
  R04<-L04[(L04[,1]==anio),]
  R05<-L05[(L05[,1]==anio),]
  R06<-L06[(L06[,1]==anio),]
  R07<-L07[(L07[,1]==anio),]
  R08<-L08[(L08[,1]==anio),]
  R09<-L09[(L09[,1]==anio),]
  R10<-L10[(L10[,1]==anio),]
  R11<-L11[(L11[,1]==anio),]
  R12<-L12[(L12[,1]==anio),]
  R13<-L13[(L13[,1]==anio),]
  R14<-L14[(L14[,1]==anio),]
  R15<-L15[(L15[,1]==anio),]
  R16<-L16[(L16[,1]==anio),]
  R17<-L17[(L17[,1]==anio),]
  R18<-L18[(L18[,1]==anio),]
  R19<-L19[(L19[,1]==anio),]
  R20<-L20[(L20[,1]==anio),]
  # Juntando estaciones
  RR<-rbind(R01,R02,R03,R04,R05,R06,R07,R08,R09,R10,
            R11,R12,R13,R14,R15,R16,R17,R18,R19,R20)
  # Precipitacion ausente 0 negativo poner cero de precipitacion
  z<-is.na(RR[,4])
  RR[z,4]<-0
  x<-RR[,4]<0
  RR[x,4]<-0
  # Acumular la precipitacion por meses
  TR<-agricolae::tapply.stat(RR[,4],RR[,c(7)],sum)
  piura <- read.csv("Piura.csv",sep = ",",header = T)
  piura$precipitacion<-TR[,2]
  return(piura)
}

1.3 Librerias necesarias

library(raster)
library(sf)
library(sp)
library(spdep)
library(gstat)
library(rgdal)
library(reshape2)
library(ggplot2)
library(RColorBrewer)

2. Interpolación

1982 - Fenómeno del Niño

AÑO=1982
region=data_anio(AÑO)
ys=region$latitud;xs=region$longitud;estaciones=region$localidad # Se guardan las latitudes, longitudes y nombres de estaciones en una variable
coordinates(region)<- ~longitud + latitud
head(region)
  orden       localidad altitud precipitacion
1     1         Ayabaca    2648       1261.21
2     2          Bernal      11          2.29
3     3      Chulucanas      89        118.61
4     4          Chusis       8          5.90
5     5 Hacienda Bigote     198        298.12
6     6 Hacienda Sumaya    1975        592.30
#DEPARTAMENTOS.shp contiene las siluetas de cada departamento del Perú
peru <- st_read("DEPARTAMENTOS.shp")[0] #st_read permite la lectura de archivos en formato shapefile (shp)
Reading layer `DEPARTAMENTOS' from data source 
  `C:\Users\Alejandro\Desktop\GED\20190310-Salvador-Piura\DEPARTAMENTOS.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 25 features and 79 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -81.32823 ymin: -18.35093 xmax: -68.65228 ymax: -0.03860597
CRS:           NA
piura<-peru[20,1] #Silueta de Piura
plot(piura) #Grafica de la silueta de piura
axis(1) #Eje x
axis(2) #Eje y

e<-extent(piura) # se toman los valores minimo y maximo de las coordenadas
grd<-expand.grid(x=seq(from=e[1],to=e[2],by=0.01), y=seq(from=e[3],to=e[4],by=0.01)) #expand.grid crea un dataframe de puntos (malla rectangular), se ingresa como parametros los valores extremos ya calculados con un intervalo entre punto y punto de 0.01
coordinates(grd)<- ~x + y #Formato de coordenadas 
gridded(grd)<- TRUE # Explica que los datos espaciales estan cuadriculados
crs(grd)<-crs(shapefile)
nuevo<-idw(precipitacion/12 ~ 1, region,grd) # se generan los valores en funcion de la precipitación para cada punto creado
[inverse distance weighted interpolation]
nuevo<-raster(nuevo) # se crea objetos de tipo raster que representan información geoespacial en forma de una cuadrícula regular.
nuevo1<-raster::mask(nuevo,piura) #mask crea una mascara en función a las coordenadas de la malla nueva y la silueta
plot(nuevo1) # Grafico de interpolación de precipitación de la región Piura
contour(nuevo1,add=TRUE,axes=FALSE)  # se dibuja el contorno según la precipitación 
text(xs,ys,estaciones,col="blue",cex=0.5) # Se grafican las estaciones de acuerdo a sus posiciones (longitud, latitud)

image(nuevo1) # Grafico de interpolación de precipitación de la región Piura usando image
contour(nuevo1,add=TRUE,axes=FALSE) # se dibuja el contorno según la precipitación 
text(xs,ys,estaciones,col="blue",cex=0.5)  # Se grafican las estaciones de acuerdo a sus posiciones (longitud, latitud)

# Imagen generada con ggplot
ppt.df<-as.data.frame(nuevo1,xy=TRUE) %>% 
  melt(id.vars=c("x","y"))
# color en gradiente
ggplot()+
  geom_raster(data=ppt.df, aes(x=x,y=y, fill=value))+
  facet_wrap(~variable)+
  labs(x="longitud",y="latitud")+
  scale_fill_viridis_c(option = "B",   direction = -1,
                       breaks = c(20,  40,  60,  80, 100),
                       labels = c(" 20 mm", " 40 mm", " 60 mm", " 80 mm", " 100 mm"),
                       name = "Precipitación\nMensual, Año 1982")+
  ggtitle("Precipitación en Piura-Peru") + 
  geom_text(aes(x = xs, y = ys, label = estaciones), size=2.5, color = "blue") +
  theme_classic()

1983 - Fenómeno del Niño

AÑO=1983
region=data_anio(AÑO)
ys=region$latitud;xs=region$longitud;estaciones=region$localidad
coordinates(region)<- ~longitud + latitud
head(region)
  orden       localidad altitud precipitacion
1     1         Ayabaca    2648       2462.24
2     2          Bernal      11        174.60
3     3      Chulucanas      89       4137.90
4     4          Chusis       8        604.00
5     5 Hacienda Bigote     198       1753.02
6     6 Hacienda Sumaya    1975        998.50
peru <- st_read("DEPARTAMENTOS.shp")[0]
Reading layer `DEPARTAMENTOS' from data source 
  `C:\Users\Alejandro\Desktop\GED\20190310-Salvador-Piura\DEPARTAMENTOS.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 25 features and 79 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -81.32823 ymin: -18.35093 xmax: -68.65228 ymax: -0.03860597
CRS:           NA
piura<-peru[20,1]
plot(piura)
axis(1)
axis(2)

e<-extent(piura)
grd<-expand.grid(x=seq(from=e[1],to=e[2],by=0.01), y=seq(from=e[3],to=e[4],by=0.01)) # by es la distancia
coordinates(grd)<- ~x + y
gridded(grd)<- TRUE
crs(grd)<-crs(shapefile)
nuevo<-idw(precipitacion/12 ~ 1, region,grd)
[inverse distance weighted interpolation]
nuevo<-raster(nuevo)
nuevo1<-raster::mask(nuevo,piura)
plot(nuevo1)
contour(nuevo1,add=TRUE,axes=FALSE)
text(xs,ys,estaciones,col="blue",cex=0.5)

image(nuevo1)
contour(nuevo1,add=TRUE,axes=FALSE)
text(xs,ys,estaciones,col="blue",cex=0.5)

# Imagen generada con ggplot
ppt.df<-as.data.frame(nuevo1,xy=TRUE) %>% 
  melt(id.vars=c("x","y"))
# color en gradiente
ggplot()+
  geom_raster(data=ppt.df, aes(x=x,y=y, fill=value))+
  facet_wrap(~variable)+
  labs(x="longitud",y="latitud")+
  scale_fill_viridis_c(option = "B",   direction = -1,
                       breaks = c(50,  100,  150,  200, 250, 300),
                       labels = c(" 50 mm", " 100 mm", " 150 mm", " 200 mm", " 250 mm", " 300 mm"),
                       name = "Precipitación\nMensual, Año 1983")+
  ggtitle("Precipitación en Piura-Peru") + 
  geom_text(aes(x = xs, y = ys, label = estaciones), size=2.5, color = "blue") +
  theme_classic()

1990 - Sequía

AÑO=1990
region=data_anio(AÑO)
ys=region$latitud;xs=region$longitud;estaciones=region$localidad
coordinates(region)<- ~longitud + latitud
head(region)
  orden       localidad altitud precipitacion
1     1         Ayabaca    2648       1087.22
2     2          Bernal      11          4.70
3     3      Chulucanas      89         41.00
4     4          Chusis       8          1.20
5     5 Hacienda Bigote     198         60.29
6     6 Hacienda Sumaya    1975        863.10
peru <- st_read("DEPARTAMENTOS.shp")[0]
Reading layer `DEPARTAMENTOS' from data source 
  `C:\Users\Alejandro\Desktop\GED\20190310-Salvador-Piura\DEPARTAMENTOS.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 25 features and 79 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -81.32823 ymin: -18.35093 xmax: -68.65228 ymax: -0.03860597
CRS:           NA
piura<-peru[20,1]
plot(piura)
axis(1)
axis(2)

e<-extent(piura)
grd<-expand.grid(x=seq(from=e[1],to=e[2],by=0.01), y=seq(from=e[3],to=e[4],by=0.01)) # by es la distancia
coordinates(grd)<- ~x + y
gridded(grd)<- TRUE
crs(grd)<-crs(shapefile)
nuevo<-idw(precipitacion/12 ~ 1, region,grd)
[inverse distance weighted interpolation]
nuevo<-raster(nuevo)
nuevo1<-raster::mask(nuevo,piura)
plot(nuevo1)
contour(nuevo1,add=TRUE,axes=FALSE)
text(xs,ys,estaciones,col="blue",cex=0.5)

image(nuevo1)
contour(nuevo1,add=TRUE,axes=FALSE)
text(xs,ys,estaciones,col="blue",cex=0.5)

# Imagen generada con ggplot
ppt.df<-as.data.frame(nuevo1,xy=TRUE) %>% 
  melt(id.vars=c("x","y"))
# color en gradiente
ggplot()+
  geom_raster(data=ppt.df, aes(x=x,y=y, fill=value))+
  facet_wrap(~variable)+
  labs(x="longitud",y="latitud")+
  scale_fill_viridis_c(option = "B",   direction = -1,
                       breaks = c(20,40,60,80),
                       labels = c(" 20 mm", " 40 mm", " 60 mm", " 80 mm"),
                       name = "Precipitación\nMensual, Año 1990")+
  ggtitle("Precipitación en Piura-Peru") + 
  geom_text(aes(x = xs, y = ys, label = estaciones), size=2.5, color = "blue") +
  theme_classic()

1991 - Sequía

AÑO=1991
region=data_anio(AÑO)
ys=region$latitud;xs=region$longitud;estaciones=region$localidad
coordinates(region)<- ~longitud + latitud
head(region)
  orden       localidad altitud precipitacion
1     1         Ayabaca    2648        977.80
2     2          Bernal      11         13.30
3     3      Chulucanas      89          0.00
4     4          Chusis       8          0.80
5     5 Hacienda Bigote     198        111.51
6     6 Hacienda Sumaya    1975        489.00
peru <- st_read("DEPARTAMENTOS.shp")[0]
Reading layer `DEPARTAMENTOS' from data source 
  `C:\Users\Alejandro\Desktop\GED\20190310-Salvador-Piura\DEPARTAMENTOS.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 25 features and 79 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -81.32823 ymin: -18.35093 xmax: -68.65228 ymax: -0.03860597
CRS:           NA
piura<-peru[20,1]
plot(piura)
axis(1)
axis(2)

e<-extent(piura)
grd<-expand.grid(x=seq(from=e[1],to=e[2],by=0.01), y=seq(from=e[3],to=e[4],by=0.01)) # by es la distancia
coordinates(grd)<- ~x + y
gridded(grd)<- TRUE
crs(grd)<-crs(shapefile)
nuevo<-idw(precipitacion/12 ~ 1, region,grd)
[inverse distance weighted interpolation]
nuevo<-raster(nuevo)
nuevo1<-raster::mask(nuevo,piura)
plot(nuevo1)
contour(nuevo1,add=TRUE,axes=FALSE)
text(xs,ys,estaciones,col="blue",cex=0.5)

image(nuevo1)
contour(nuevo1,add=TRUE,axes=FALSE)
text(xs,ys,estaciones,col="blue",cex=0.5)

# Imagen generada con ggplot
ppt.df<-as.data.frame(nuevo1,xy=TRUE) %>% 
  melt(id.vars=c("x","y"))
# color en gradiente
ggplot()+
  geom_raster(data=ppt.df, aes(x=x,y=y, fill=value))+
  facet_wrap(~variable)+
  labs(x="longitud",y="latitud")+
  scale_fill_viridis_c(option = "B",   direction = -1,
                       breaks = c(20,40,60,80),
                       labels = c(" 20 mm", " 40 mm", " 60 mm", " 80 mm"),
                       name = "Precipitación\nMensual, Año 1991")+
  ggtitle("Precipitación en Piura-Peru") + 
  geom_text(aes(x = xs, y = ys, label = estaciones), size=2.5, color = "blue") +
  theme_classic()

1993 - Sequía

AÑO=1993
region=data_anio(AÑO)
ys=region$latitud;xs=region$longitud;estaciones=region$localidad
coordinates(region)<- ~longitud + latitud
head(region)
  orden       localidad altitud precipitacion
1     1         Ayabaca    2648       1473.90
2     2          Bernal      11         22.10
3     3      Chulucanas      89          0.00
4     4          Chusis       8          0.00
5     5 Hacienda Bigote     198        786.13
6     6 Hacienda Sumaya    1975       1206.20
peru <- st_read("DEPARTAMENTOS.shp")[0]
Reading layer `DEPARTAMENTOS' from data source 
  `C:\Users\Alejandro\Desktop\GED\20190310-Salvador-Piura\DEPARTAMENTOS.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 25 features and 79 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -81.32823 ymin: -18.35093 xmax: -68.65228 ymax: -0.03860597
CRS:           NA
piura<-peru[20,1]
plot(piura)
axis(1)
axis(2)

e<-extent(piura)
grd<-expand.grid(x=seq(from=e[1],to=e[2],by=0.01), y=seq(from=e[3],to=e[4],by=0.01)) # by es la distancia
coordinates(grd)<- ~x + y
gridded(grd)<- TRUE
crs(grd)<-crs(shapefile)
nuevo<-idw(precipitacion/12 ~ 1, region,grd)
[inverse distance weighted interpolation]
nuevo<-raster(nuevo)
nuevo1<-raster::mask(nuevo,piura)
plot(nuevo1)
contour(nuevo1,add=TRUE,axes=FALSE)
text(xs,ys,estaciones,col="blue",cex=0.5)

image(nuevo1)
contour(nuevo1,add=TRUE,axes=FALSE)
text(xs,ys,estaciones,col="blue",cex=0.5)

# Imagen generada con ggplot
ppt.df<-as.data.frame(nuevo1,xy=TRUE) %>% 
  melt(id.vars=c("x","y"))
# color en gradiente
ggplot()+
  geom_raster(data=ppt.df, aes(x=x,y=y, fill=value))+
  facet_wrap(~variable)+
  labs(x="longitud",y="latitud")+
  scale_fill_viridis_c(option = "B",   direction = -1,
                       breaks = c(50,100,150),
                       labels = c(" 50 mm", " 100 mm", " 150 mm"),
                       name = "Precipitación\nMensual, Año 1993")+
  ggtitle("Precipitación en Piura-Peru") + 
  geom_text(aes(x = xs, y = ys, label = estaciones), size=2.5, color = "blue") +
  theme_classic()

1994 - Sequía

AÑO=1994
region=data_anio(AÑO)
ys=region$latitud;xs=region$longitud;estaciones=region$localidad
coordinates(region)<- ~longitud + latitud
head(region)
  orden       localidad altitud precipitacion
1     1         Ayabaca    2648        216.00
2     2          Bernal      11         34.81
3     3      Chulucanas      89          0.00
4     4          Chusis       8         14.30
5     5 Hacienda Bigote     198        701.70
6     6 Hacienda Sumaya    1975       1340.80
peru <- st_read("DEPARTAMENTOS.shp")[0]
Reading layer `DEPARTAMENTOS' from data source 
  `C:\Users\Alejandro\Desktop\GED\20190310-Salvador-Piura\DEPARTAMENTOS.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 25 features and 79 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -81.32823 ymin: -18.35093 xmax: -68.65228 ymax: -0.03860597
CRS:           NA
piura<-peru[20,1]
plot(piura)
axis(1)
axis(2)

e<-extent(piura)
grd<-expand.grid(x=seq(from=e[1],to=e[2],by=0.01), y=seq(from=e[3],to=e[4],by=0.01)) # by es la distancia
coordinates(grd)<- ~x + y
gridded(grd)<- TRUE
crs(grd)<-crs(shapefile)
nuevo<-idw(precipitacion/12 ~ 1, region,grd)
[inverse distance weighted interpolation]
nuevo<-raster(nuevo)
nuevo1<-raster::mask(nuevo,piura)
plot(nuevo1)
contour(nuevo1,add=TRUE,axes=FALSE)
text(xs,ys,estaciones,col="blue",cex=0.5)

image(nuevo1)
contour(nuevo1,add=TRUE,axes=FALSE)
text(xs,ys,estaciones,col="blue",cex=0.5)

# Imagen generada con ggplot
ppt.df<-as.data.frame(nuevo1,xy=TRUE) %>% 
  melt(id.vars=c("x","y"))
# color en gradiente
ggplot()+
  geom_raster(data=ppt.df, aes(x=x,y=y, fill=value))+
  facet_wrap(~variable)+
  labs(x="longitud",y="latitud")+
  scale_fill_viridis_c(option = "B",   direction = -1,
                       breaks = c(20,40,60,80,100,120),
                       labels = c(" 20 mm", " 40 mm", " 60 mm", " 80 mm", " 100 mm", " 120 mm"),
                       name = "Precipitación\nMensual, Año 1994")+
  ggtitle("Precipitación en Piura-Peru") + 
  geom_text(aes(x = xs, y = ys, label = estaciones), size=2.5, color = "blue") +
  theme_classic()

1995- Sequía

AÑO=1995
region=data_anio(AÑO)
ys=region$latitud;xs=region$longitud;estaciones=region$localidad
coordinates(region)<- ~longitud + latitud
head(region)
  orden       localidad altitud precipitacion
1     1         Ayabaca    2648        959.80
2     2          Bernal      11         13.85
3     3      Chulucanas      89          0.00
4     4          Chusis       8         10.90
5     5 Hacienda Bigote     198        247.14
6     6 Hacienda Sumaya    1975        710.00
peru <- st_read("DEPARTAMENTOS.shp")[0]
Reading layer `DEPARTAMENTOS' from data source 
  `C:\Users\Alejandro\Desktop\GED\20190310-Salvador-Piura\DEPARTAMENTOS.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 25 features and 79 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -81.32823 ymin: -18.35093 xmax: -68.65228 ymax: -0.03860597
CRS:           NA
piura<-peru[20,1]
plot(piura)
axis(1)
axis(2)

e<-extent(piura)
grd<-expand.grid(x=seq(from=e[1],to=e[2],by=0.01), y=seq(from=e[3],to=e[4],by=0.01)) # by es la distancia
coordinates(grd)<- ~x + y
gridded(grd)<- TRUE
crs(grd)<-crs(shapefile)
nuevo<-idw(precipitacion/12 ~ 1, region,grd)
[inverse distance weighted interpolation]
nuevo<-raster(nuevo)
nuevo1<-raster::mask(nuevo,piura)
plot(nuevo1)
contour(nuevo1,add=TRUE,axes=FALSE)
text(xs,ys,estaciones,col="blue",cex=0.5)

image(nuevo1)
contour(nuevo1,add=TRUE,axes=FALSE)
text(xs,ys,estaciones,col="blue",cex=0.5)

# Imagen generada con ggplot
ppt.df<-as.data.frame(nuevo1,xy=TRUE) %>% 
  melt(id.vars=c("x","y"))
# color en gradiente
ggplot()+
  geom_raster(data=ppt.df, aes(x=x,y=y, fill=value))+
  facet_wrap(~variable)+
  labs(x="longitud",y="latitud")+
  scale_fill_viridis_c(option = "B",   direction = -1,
                       breaks = c(20,40,60),
                       labels = c(" 20 mm", " 40 mm", " 60 mm"),
                       name = "Precipitación\nMensual, Año 1995")+
  ggtitle("Precipitación en Piura-Peru") + 
  geom_text(aes(x = xs, y = ys, label = estaciones), size=2.5, color = "blue") +
  theme_classic()

1996 - Sequía

AÑO=1996
region=data_anio(AÑO)
ys=region$latitud;xs=region$longitud;estaciones=region$localidad
coordinates(region)<- ~longitud + latitud
head(region)
  orden       localidad altitud precipitacion
1     1         Ayabaca    2648        816.70
2     2          Bernal      11          5.53
3     3      Chulucanas      89          0.00
4     4          Chusis       8          3.80
5     5 Hacienda Bigote     198         77.58
6     6 Hacienda Sumaya    1975        753.60
peru <- st_read("DEPARTAMENTOS.shp")[0]
Reading layer `DEPARTAMENTOS' from data source 
  `C:\Users\Alejandro\Desktop\GED\20190310-Salvador-Piura\DEPARTAMENTOS.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 25 features and 79 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -81.32823 ymin: -18.35093 xmax: -68.65228 ymax: -0.03860597
CRS:           NA
piura<-peru[20,1]
plot(piura)
axis(1)
axis(2)

e<-extent(piura)
grd<-expand.grid(x=seq(from=e[1],to=e[2],by=0.01), y=seq(from=e[3],to=e[4],by=0.01)) # by es la distancia
coordinates(grd)<- ~x + y
gridded(grd)<- TRUE
crs(grd)<-crs(shapefile)
nuevo<-idw(precipitacion/12 ~ 1, region,grd)
[inverse distance weighted interpolation]
nuevo<-raster(nuevo)
nuevo1<-raster::mask(nuevo,piura)
plot(nuevo1)
contour(nuevo1,add=TRUE,axes=FALSE)
text(xs,ys,estaciones,col="blue",cex=0.5)

image(nuevo1)
contour(nuevo1,add=TRUE,axes=FALSE)
text(xs,ys,estaciones,col="blue",cex=0.5)

# Imagen generada con ggplot
ppt.df<-as.data.frame(nuevo1,xy=TRUE) %>% 
  melt(id.vars=c("x","y"))
# color en gradiente
ggplot()+
  geom_raster(data=ppt.df, aes(x=x,y=y, fill=value))+
  facet_wrap(~variable)+
  labs(x="longitud",y="latitud")+
  scale_fill_viridis_c(option = "B",   direction = -1,
                       breaks = c(10,20,30,40,50,60),
                       labels = c(" 10 mm"," 20 mm"," 30 mm", " 40 mm"," 50 mm", " 60 mm"),
                       name = "Precipitación\nMensual, Año 1996")+
  ggtitle("Precipitación en Piura-Peru") + 
  geom_text(aes(x = xs, y = ys, label = estaciones), size=2.5, color = "blue") +
  theme_classic()

1997 - Fenómeno del Niño

AÑO=1997
region=data_anio(AÑO)
ys=region$latitud;xs=region$longitud;estaciones=region$localidad
coordinates(region)<- ~longitud + latitud
head(region)
  orden       localidad altitud precipitacion
1     1         Ayabaca    2648       1223.00
2     2          Bernal      11         54.63
3     3      Chulucanas      89        142.33
4     4          Chusis       8         55.70
5     5 Hacienda Bigote     198        407.50
6     6 Hacienda Sumaya    1975        789.70
peru <- st_read("DEPARTAMENTOS.shp")[0]
Reading layer `DEPARTAMENTOS' from data source 
  `C:\Users\Alejandro\Desktop\GED\20190310-Salvador-Piura\DEPARTAMENTOS.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 25 features and 79 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -81.32823 ymin: -18.35093 xmax: -68.65228 ymax: -0.03860597
CRS:           NA
piura<-peru[20,1]
plot(piura)
axis(1)
axis(2)

e<-extent(piura)
grd<-expand.grid(x=seq(from=e[1],to=e[2],by=0.01), y=seq(from=e[3],to=e[4],by=0.01)) # by es la distancia
coordinates(grd)<- ~x + y
gridded(grd)<- TRUE
crs(grd)<-crs(shapefile)
nuevo<-idw(precipitacion/12 ~ 1, region,grd)
[inverse distance weighted interpolation]
nuevo<-raster(nuevo)
nuevo1<-raster::mask(nuevo,piura)
plot(nuevo1)
contour(nuevo1,add=TRUE,axes=FALSE)
text(xs,ys,estaciones,col="blue",cex=0.5)

image(nuevo1)
contour(nuevo1,add=TRUE,axes=FALSE)
text(xs,ys,estaciones,col="blue",cex=0.5)

# Imagen generada con ggplot
ppt.df<-as.data.frame(nuevo1,xy=TRUE) %>% 
  melt(id.vars=c("x","y"))
# color en gradiente
ggplot()+
  geom_raster(data=ppt.df, aes(x=x,y=y, fill=value))+
  facet_wrap(~variable)+
  labs(x="longitud",y="latitud")+
  scale_fill_viridis_c(option = "B",   direction = -1,
                       breaks = c(20,  40,  60,  80, 100),
                       labels = c(" 20 mm", " 40 mm", " 60 mm", " 80 mm", " 100 mm"),
                       name = "Precipitación\nMensual, Año 1997")+
  ggtitle("Precipitación en Piura-Peru") + 
  geom_text(aes(x = xs, y = ys, label = estaciones), size=2.5, color = "blue") +
  theme_classic()

1998 - Fenómeno del Niño

AÑO=1998
region=data_anio(AÑO)
ys=region$latitud;xs=region$longitud;estaciones=region$localidad
coordinates(region)<- ~longitud + latitud
head(region)
  orden       localidad altitud precipitacion
1     1         Ayabaca    2648       1990.00
2     2          Bernal      11       1200.64
3     3      Chulucanas      89        333.50
4     4          Chusis       8        983.10
5     5 Hacienda Bigote     198       2107.40
6     6 Hacienda Sumaya    1975       1474.10
peru <- st_read("DEPARTAMENTOS.shp")[0]
Reading layer `DEPARTAMENTOS' from data source 
  `C:\Users\Alejandro\Desktop\GED\20190310-Salvador-Piura\DEPARTAMENTOS.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 25 features and 79 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -81.32823 ymin: -18.35093 xmax: -68.65228 ymax: -0.03860597
CRS:           NA
piura<-peru[20,1]
plot(piura)
axis(1)
axis(2)

e<-extent(piura)
grd<-expand.grid(x=seq(from=e[1],to=e[2],by=0.01), y=seq(from=e[3],to=e[4],by=0.01)) # by es la distancia
coordinates(grd)<- ~x + y
gridded(grd)<- TRUE
crs(grd)<-crs(shapefile)
nuevo<-idw(precipitacion/12 ~ 1, region,grd)
[inverse distance weighted interpolation]
nuevo<-raster(nuevo)
nuevo1<-raster::mask(nuevo,piura)
plot(nuevo1)
contour(nuevo1,add=TRUE,axes=FALSE)
text(xs,ys,estaciones,col="blue",cex=0.5)

image(nuevo1)
contour(nuevo1,add=TRUE,axes=FALSE)
text(xs,ys,estaciones,col="blue",cex=0.5)

# Imagen generada con ggplot
ppt.df<-as.data.frame(nuevo1,xy=TRUE) %>% 
  melt(id.vars=c("x","y"))
# color en gradiente
ggplot()+
  geom_raster(data=ppt.df, aes(x=x,y=y, fill=value))+
  facet_wrap(~variable)+
  labs(x="longitud",y="latitud")+
  scale_fill_viridis_c(option = "B",   direction = -1,
                       breaks = c(50,  100,  150,  200, 250),
                       labels = c(" 50 mm", " 100 mm", " 150 mm", " 200 mm", " 250 mm"),
                       name = "Precipitación\nMensual, Año 1998")+
  ggtitle("Precipitación en Piura-Peru") + 
  geom_text(aes(x = xs, y = ys, label = estaciones), size=2.5, color = "blue") +
  theme_classic()

1999 - Fenómeno de la Niña

AÑO=1999
region=data_anio(AÑO)
ys=region$latitud;xs=region$longitud;estaciones=region$localidad
coordinates(region)<- ~longitud + latitud
head(region)
  orden       localidad altitud precipitacion
1     1         Ayabaca    2648       2103.00
2     2          Bernal      11         71.26
3     3      Chulucanas      89        431.23
4     4          Chusis       8         48.20
5     5 Hacienda Bigote     198        845.10
6     6 Hacienda Sumaya    1975       1862.50
peru <- st_read("DEPARTAMENTOS.shp")[0]
Reading layer `DEPARTAMENTOS' from data source 
  `C:\Users\Alejandro\Desktop\GED\20190310-Salvador-Piura\DEPARTAMENTOS.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 25 features and 79 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -81.32823 ymin: -18.35093 xmax: -68.65228 ymax: -0.03860597
CRS:           NA
piura<-peru[20,1]
plot(piura)
axis(1)
axis(2)

e<-extent(piura)
grd<-expand.grid(x=seq(from=e[1],to=e[2],by=0.01), y=seq(from=e[3],to=e[4],by=0.01)) # by es la distancia
coordinates(grd)<- ~x + y
gridded(grd)<- TRUE
crs(grd)<-crs(shapefile)
nuevo<-idw(precipitacion/12 ~ 1, region,grd)
[inverse distance weighted interpolation]
nuevo<-raster(nuevo)
nuevo1<-raster::mask(nuevo,piura)
plot(nuevo1)
contour(nuevo1,add=TRUE,axes=FALSE)
text(xs,ys,estaciones,col="blue",cex=0.5)

image(nuevo1)
contour(nuevo1,add=TRUE,axes=FALSE)
text(xs,ys,estaciones,col="blue",cex=0.5)

# Imagen generada con ggplot
ppt.df<-as.data.frame(nuevo1,xy=TRUE) %>% 
  melt(id.vars=c("x","y"))
# color en gradiente
ggplot()+
  geom_raster(data=ppt.df, aes(x=x,y=y, fill=value))+
  facet_wrap(~variable)+
  labs(x="longitud",y="latitud")+
  scale_fill_viridis_c(option = "B",   direction = -1,
                       breaks = c(50,  100,  150),
                       labels = c(" 50 mm", " 100 mm", " 150 mm"),
                       name = "Precipitación\nMensual, Año 1999")+
  ggtitle("Precipitación en Piura-Peru") + 
  geom_text(aes(x = xs, y = ys, label = estaciones), size=2.5, color = "blue") +
  theme_classic()

2000 - Fenómeno de la Niña

AÑO=2000
region=data_anio(AÑO)
ys=region$latitud;xs=region$longitud;estaciones=region$localidad
coordinates(region)<- ~longitud + latitud
head(region)
  orden       localidad altitud precipitacion
1     1         Ayabaca    2648       1599.43
2     2          Bernal      11         16.09
3     3      Chulucanas      89        287.91
4     4          Chusis       8         26.70
5     5 Hacienda Bigote     198        818.40
6     6 Hacienda Sumaya    1975       1599.20
peru <- st_read("DEPARTAMENTOS.shp")[0]
Reading layer `DEPARTAMENTOS' from data source 
  `C:\Users\Alejandro\Desktop\GED\20190310-Salvador-Piura\DEPARTAMENTOS.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 25 features and 79 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -81.32823 ymin: -18.35093 xmax: -68.65228 ymax: -0.03860597
CRS:           NA
piura<-peru[20,1]
plot(piura)
axis(1)
axis(2)

e<-extent(piura)
grd<-expand.grid(x=seq(from=e[1],to=e[2],by=0.01), y=seq(from=e[3],to=e[4],by=0.01)) # by es la distancia
coordinates(grd)<- ~x + y
gridded(grd)<- TRUE
crs(grd)<-crs(shapefile)
nuevo<-idw(precipitacion/12 ~ 1, region,grd)
[inverse distance weighted interpolation]
nuevo<-raster(nuevo)
nuevo1<-raster::mask(nuevo,piura)
plot(nuevo1)
contour(nuevo1,add=TRUE,axes=FALSE)
text(xs,ys,estaciones,col="blue",cex=0.5)

image(nuevo1)
contour(nuevo1,add=TRUE,axes=FALSE)
text(xs,ys,estaciones,col="blue",cex=0.5)

# Imagen generada con ggplot
ppt.df<-as.data.frame(nuevo1,xy=TRUE) %>% 
  melt(id.vars=c("x","y"))
# color en gradiente
ggplot()+
  geom_raster(data=ppt.df, aes(x=x,y=y, fill=value))+
  facet_wrap(~variable)+
  labs(x="longitud",y="latitud")+
  scale_fill_viridis_c(option = "B",   direction = -1,
                       breaks = c(20,  40,  60,  80, 100, 120),
                       labels = c(" 20 mm", " 40 mm", " 60 mm", " 80 mm", " 100 mm", " 120 mm"),
                       name = "Precipitación\nMensual, Año 2000")+
  ggtitle("Precipitación en Piura-Peru") + 
  geom_text(aes(x = xs, y = ys, label = estaciones), size=2.5, color = "blue") +
  theme_classic()